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ABSTRACT 



Context. Outflows of young stars drive shocks into dusty, molecular regions. Most models of such shocks are restricted by the 
assumptions that they are steady and propagating in directions perpendicular to the magnetic fields. However, the media through 
which shocks propagate are inhomogeneous and shocks are not steady. Furthermore, only a small fraction of shocks are nearly 
perpendicular. 

Aims. We identify features that develop when a shock encounters a density inhomogeneity and ascertain if any part of the precursor 
region of a non-steady multifluid shock ever behaves in a quasi-steady fashion. If it does, some time-dependent shocks may be 
modelled approximately without solving the time-dependent hydromagnetic equations. 

Methods. We use the code employed previously to produce the first time-dependent simulations of fast-mode oblique C-type shocks 
including a self-consistent calculation of the thermal and ionisation balances and a fluid treatment of grains. 

Results. Simulations were made for initially steady oblique C-type shocks, each of which encounters one of three types of density 
inhomogeneities. For a semi-finite inhomogeneity with a density larger than the surrounding medium's, a transmitted shock evolves 
from being of J-type to a steady C-type shock on a timescale comparable to the ion-flow time through it. A sufficiently upstream 
part of the precursor of an evolving J-type shock is quasi- steady. The ion-flow timescale is also relevant for the evolution of a shock 
moving into a region of decreasing density. The models for shocks propagating into regions in which the density increases and then 
decreases to its initial value cannot be entirely described in terms of the results obtained for monotonically increasing and decreasing 
densities. 

Conclusions. We present the first time-dependent simulations of dusty C-type shocks interacting with density perturbations. We 
studied the transient evolution of the shock structure and find that the initial interaction always produces a transition to a J-type shock. 
Furthermore, the long-term evolution back to a C-type shock cannot always be approximated by quasi- steady models. 

Key words. MHD - shock waves - ISM: dust - ISM: jets and outflows 



1. Introduction 

Jets and winds associated with low and high mass proto- 
stellar objects interact w i th the molecular clouds surround- 
ing them (e.g. lArce et al" 1 1200% . Shocks are driven into the 
clouds sweeping up cloud material and producing large scale 
(0.1 - 1 pc) molecular outflows wit h velocities of 10- 1 00 km 
s" 1 (e.g. iBachiller & TafallaL 1 19991: iReipurth & BalR 120011: 
ISantiago-Garcia etal ■L 120091) . These outflows are often bipo- 
lar and are extreme ly common in young, low- mass stars 
(iBallv & LadaL fl983l) .~ as well as high-mass star s (ICesaroniL 
120051: IShepherdL 120051: lL6pez-Sepulcre et alll2009h . 

A s the fractional io nisation in molecular clouds is low (x < 
10" 6 ; iDalgarnol 120061 and reference therein), the neutral gas 
and magnetic field are weakly coupled. This coupling is me- 
diated via ion-neutral and grain-neutral collisions, as the mag- 
netic field forces the charged particles to move through the bulk 
of neutral particles. When a shock forms, the ion-neutral colli- 
sions accelerate, compress and heat the neutral gas upstream of 
the shock, forming a magnetic precursor. If the shock velocity 
is low and the postshock cooling efficient, the shock becomes 
continuous in all fluids and is then referred to as a C-type shock 
dMullanl [l97lt iDrainei fl980l) . Observational studies of the ve- 



locity widths of SiO emission suggest that most shocks near low 
mass proto- stellar objects are C-type (e.g. lMartin-Pintado et all 
Il992h . 

Previous studies of C-type shocks in dusty plasmas have 
focu sed on the shock structure under steady state conditions 
(e.g. iDraine et all 119831: iPffipp et all Il990l: iPilipp & HartquistL 
119941: IWardlel 119981: iGuillet et all 12007b . However, SiO obser- 
vations suggest that proto- stellar jets and wind s interact with 
clumpy structures a l ong their propagation axis ( Mikami et all 
119921: iLefloch et all 119981: iJorgensen et all |2004l). Therefore 
the ste ady state assumption must be relaxed. ICiolek & Robergd 
(120021) were the firsQ to use a time-dependent multifluid mag- 
netohydrodynamics (MHD) code, including dust grains dynam- 
ics, to study fast-mode C-type shocks. However, their study is 
limited to perpendicu lar shocks. The multifluid approach devel- 
oped by Falle ( 2003) ov ercom es this restriction. Using the Falle 
method, IVan Loo et al.l (I2Q09L hereafter Paper I) produced the 
first time-dependent simulations of fast-mode, oblique C-type 
shocks that evolve to steady state and include a selfconsistent 
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1 Time-depende nt multifluid shock stu dies prior to 

ICiolek & Robergd 120021) , e.g. those oflTothldl994n :I Smith & Mac Low! 
d!997l) : IStond dl997l) and IChieze et all d!998l) . did not include grain 
dynamics. 
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calculation of the thermal and ionisation balances and a fluid 
treatment of grain dynamics. 

While the time-dependent simulations of Paper I still focus 
on steady state shock structures, we now extend this work by 
studying the temporal evolution of fast-mode C-type shocks in- 
teracting with regions of inhomogeneous density. In Sect. [2 we 
describe the numerical model with the initial conditions. We then 
apply the code to an oblique shock propagating into a region of 
varying density (Sect. [3]) and, in Sect. 01 we discuss these results 
and give our conclusions. 

2. The model 

2.1. Numerical code 

The numerical code that we use in this paper follows the dynam- 
ics of a multifluid plasma consisting of neutrals, ions, electrons 
and N dust grain species. Here we only consider a single fluid of 
spherical, large grains with a radius and mass of a g = 0.4/mi and 
8.04 x 10" 13 g. (The density of the grains is assumed to be 3 g 
cm -3 .) We have included relevant mass transfer processes, such 
as electron recombination with Mg + and dissociative recombina- 
tion with HCO + , and radiative cooling by O, CO, H 2 and H 2 0. 
Also, the mass and charge transfer from ions and electrons to the 
dust grains is taken into account to calculate the average grain 
charge. This way we calculate self-consistently the thermal and 
ionisation balances including appropriate microphysics. 

Th e details of the numerical scheme, which is based on lFallel 
(I2003I) . are given in Paper I. As in Paper I, we advance the 
scheme explicitly, e ven though t his implies a restriction on the 
stable time- step, i.e. lFallel (I2003I) 



where r a( i is the ambipolar resistivity, r H the Hall resistivity and 
Ax the cell spacing. Thus, the time step becomes small for a high 
numerical resolution, or when the Hall resistivity becomes large 
compared to the ambipolar resistivity. In our simulations the Hall 
resistivity is only ever comparable to the ambipolar resistivity 
implying that the stable time-step is not severely restricted. 

The adopted c hemistry is based on the network given in 
iPilipp et all (Il990h . and we used the rate coefficients that they 
did. The code allows the option of incorporating as many ad- 
vected scalar equations with source terms (i. e. rate equations 
for fractional abundances) as necessary and solving them with 
the same methods as those employed to solve the fluid equa- 
tions. We used rate equations and charge neutrality to obtain the 
fractional abundances of ions and electrons and the average grain 
charge. We assumed that cosmic ray induced ionisation leads im- 
mediately to the production of heavy molecular ions (e.g. HCO + 
and H30 + ) which involves the justifiable neglect of H3 + recom- 
bination. The branching ratio between HCO + and HsO + forma- 
tion depends on the abundances of simple neutral species like 
O and CO and the rate coefficients of the reaction rates of H 3 + 
with them. Charge transfer with metals (e.g. Mg) creates metallic 
ions. All gas phase ions re combine with electrons and on grains. 
Unlike lPilipp et alJ (Il990l) . we did not include the possibility that 
Mg + reacts with H 2 ; it does so only in regions with sufficiently 
high temperatures that the grain charge is high enough and the 
electron fractional abundance low enough that recombination is 
dominated by recombination with grains leading to a negligible 
difference between the rates at which Mg + and MgH + are re- 
moved. While the abundances of ions and electrons and grain 
charge can be studied easily with rate equations, the conversion 



of O to H 2 is more problematic as the relevant reaction rates in- 
crease by many orders of magnitude as shock gas is heated from 
several hundred to a thousand degrees. For an H 2 number density 
of 10 5 cm -3 , the timescale for the removal of O in reactions with 
H 2 is 1.4 xlO 12 s, 1.4 xlO 9 s, and 6.7 xlO 7 s at 300K, 600K, 
and 1000K, respectively. The corresponding removal timescale 
of OH in reactions with H 2 are 1.5 xlO 9 s, 3.2 xlO 7 s, and 4.2 
xlO 6 s. We could use the code to deal with the O, OH, and H 2 
chemistry even at temperatures in excess of 10 3 K if we placed 
upper limits on the reaction rates; this would still give oxygen 
chemistry that would be close to the appropriate high tempera- 
ture equilibrium. However, here we simply assumed that almost 
all oxygen not in CO is in O until T = 500K. After it reached 
500K almost all oxygen not in CO was taken to be in H 2 0. It 
will remain mostly in H 2 for very long times even after the 
postshock gas has cooled. This approach led to neutral temper- 
atures in the 300-500K range that are somewhat higher and in 
the 500-700K range that are somewhat lower than a more thor- 
ough treatment would have. However, this results in a negligible 
difference in the shock behaviour which is much more strongly 
influenced by the ionisation structure, which in turn is regulated 
in this regime primarily by cosmic ray induced ionization and 
recombination on grains, and grain dynamics than the details of 
the neutral gas temperature structure. 

2.2. Initial conditions 

Unlike the simulations described in Paper I, for which the initial 
conditions are J-type shocks, each of our simulations starts from 
a steady C-type shock propagating through a medium with up- 
stream hydrogen nuclei number density of either nu = 10 4 cm~ 3 
or nu = 10 5 cm~ 3 . The shock moves parallel to the x-axis with a 
velocity V s of 25 km s" 1 . For both cases the initial upstream fluid 
temperatures are 8.4K and the upstream magnetic field strengths 
are B = 10" 4 G. The magnetic field lies in the x, z/-plane with the 
angle between the shock normal and the field being 45°. 

As we are interested in the interaction of a C-type shock with 
an inhomogeneous medium, we introduce an upstream density 
perturbation. The shape of the perturbation is chosen so that it 
reasonably reproduces all or part of the densi ty profile of ob- 
served clumps or cores (e.g. iTaf alia et al .1 12004) . Each perturba- 
tion is described solely by its width and the density contrast (i.e. 
a multiple of initial upstream density). Unless otherwise stated, 
the width and density contrast are 10 17 cm and 10 respectively. 

2.3. Computational considerations 
2.3.1. Computational grid 

Some care must be given to the size of the computational grid for 
models with varying upstream density. This is because a change 
in the upstream density alters the shock width. For all our simu- 
lations, the constraint 

is met (with x me fractional ionisation). The shock thickness 
is then determined by grain-neutral rather than ion-neutral drag 
(see Paper I) and given approximately by 

t 01 m24 / V s r 1 / v A \/10- 2 cm" 3 \ 

L = 2 ' 1Xl0 ^(knT^) (kn^)(^— ) Cm ' (1) 

where va is the Alfven speed and n, the ion number density. The 
shock width in our simulations varies by an order of a mag- 
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nitude. We need to ensure that the numerical domain is large 
enough to contain the largest shock width but has enough cells 
to resolve structures on the finest scale, i.e the shock width at the 
highest density. Therefore, we estimate the shock thickness on 
both densities using Eq. [T] The numerical domain is then set to 
-10L < x < 10L with L the shock width at the lowest density. 
As the shock width at the highest density is roughly an order of 
magnitude smaller and we want to resolve it with approximately 
10 cells across the shock structure, we adopt a uniform grid of 
1200 cells for the numerical domain. We impose a fixed inflow 
at the upstream boundary and a free-flow boundary condition at 
the downstream side. 

2.3.2. Varying shock speed 

When a shock moves through a region of varying density, its 
propagation speed changes. As we follow the evolution in the 
initial shock frame, the new shock structures (see Sect. [3]) will 
eventually move off the grid. In order to follow their evolution 
properly, we adjust the velocity of the initial shock frame with a 
Lorentz transformation along the x-axis. 

2.3.3. Charged fluid inertia 

Although our initial conditions are given by a C-type shock, 
the interaction with a density inhomogeneity produces a neu- 
tral subshock (see Sect. [3]). Within the subshock, our assump- 
tion that the charged fluid inertia is negligible no longer holds. 
However, the neutral subshock struct ure can be a ccurately mod- 
elled with our numerical model (see lFalldJ2003l) . Furthermore, 
the inertial phase is expected to be shortlived compared to the 
timescales of the shock-clump interaction and the subsequent 
approach to steady state. We can expect that the evolution of the 
shock shortly after the interaction is modified by charged fluid 
inertia, but that the changes are modest. 

3. Results 

Several aspects of the shock interaction with an inhomogeneity 
are consistent for all our simulations. We discuss these first be- 
fore commenting on specific models. 

In each of our models the initial shock width is larger than 
or comparable to the length scale of the perturbation. The in- 
teraction with the perturbation then leads to a transmitted wave- 
reflected wave pair separated by a contact discontinuity. Note 
that, if the perturbation length scale were large compared to the 
shock width, there would not be a reflected wave. The relative 
strength of the transmitted and reflected components is deter- 
mined by both the initial shock speed and the amplitude of the 
perturbation. 

The transmitted shock is initially a J-type shock with a pre- 
cursor and a subshock. The upstream conditions change faster 
than information can propagate across the front of the subshock. 
After some time, collisions between ions and neutrals generate 
a neutral precursor which then evolves into a C-type shock. The 
transmitted C-type shock is narrower than the initial shock if the 
upstream density increases and broader if the upstream density 
decreases. 

While the transmitted wave is always a shock, the reflected 
wave which propagates into the post-shock flow can either be a 
rarefaction wave or a shock. As the reflected wave connects the 
far downstream flow with the downstream flow of the transmit- 
ted shock, it is clear that a rarefaction forms when the post- shock 



density of the transmitted shock is lower than the far downstream 
density and a shock when it is greater. Similarly to the transmit- 
ted shock, the reflected shock evolves from a J-type shock into a 
C-type shock. 

In this paper we focus primarily on the transmitted shocks. 
We do not follow the long term evolution of the reflected 
waves. The reflected component thus eventually moves off the 
grid at the downstream boundary. However, numerical artefacts 
arise when shock s are reflected from free-flow boundaries (e.g. 
lHedstronifl979h . Although we can use non-reflecting boundary 
conditions to overcome this problem, we prefer to simply re- 
move the reflected wave from the computational domain. This 
can be done because the reflected wave and transmitted shock 
are separated by a contact discontinuity. 

3.1. Increasing density 

We first follow the evolution of a C-type shock in the molecular 
outflow from a proto- stellar object moving into a denser region. 
We consider a 25 km s" 1 shock propagating from a region of 
Ha = 10 4 cm~ 3 into a higher density plasma. The initial shock 
is a strong fast-mode shock with an Alfvenic Mach number of 
« 11.5. 

When the C-type shock encounters the density perturbation, 
two J-type shocks, a reflected and a transmitted one, form within 
3 x 10 3 yr. The neutral subshocks can be easily seen in both the 
velocity and temperature plots of Figs. [T] [2]and[3j i.e. the dis- 
continuous jump in the neutral velocity along the x-axis and the 
drastic cooling of the different fluids. The temperature plots also 
shows that the transmitted shock is considerably stronger than 
the reflected one, as the maximum temperature in the transmit- 
ted shock is roughly 5 times higher. The temperatures of both 
the transmitted and reflected shocks decrease at later times. This 
is due to a reduced relative velocity between the neutrals and 
the charged particles as the energy transfer by collisions heats 
the gas. The lower temperatures and relative velocities suggest 
that the higher upstream density regime is less effective for grain 
sputtering. This agrees with the results o f lMavetal.1 (120001) who 
reported that the elemental fraction of Si sputtered from olivine 
increased with shock velocity, but at a given shock velocity, 
dropped as the pre-shock density increased. 

In these early stages of the evolution we can already see the 
effect of the higher density on the charged particles fluid and 
the magnetic field. In the high density region the Hall conductiv- 
ity becomes larger than the Pedersen conductivity (see definition 
in Paper I). This changes the wave propagation upstream. This 
is associated with a change of the struct ure of the evolu tionary 
trajectory in the B y - B z phase space (cf. IWardla 1 19981) . When 
the Hall conductivity exceeds the Pedersen conductivity, a spiral 
node develops in the trajectory for those regions corresponding 
to the flow around and just upstream of the leading part of the 
shock pre cursor The s hock structure then contains large oscilla- 
tions (see iFallei 120031 and Paper I) as seen in the x-component 
of the ion and electron velocities in Fig. [2] 

The J-type phase of the shocks is quite short. Within 7.92 x 
10 4 yr after impact, corresponding to 2 ion shock crossing times, 
the transmitted shock structure has become a stable C-type 
shock. For the reflected shock the time for the shock to become a 
stable C-type shock is somewhat shorter. The transmitted C-type 
shock is much narrower than the initial shock (roughly 1/8 of the 
original shock width). It also has a lower propagation speed of 
16.3 km s" 1 as the shock slows down when moving into a region 
of higher density and pressure. 
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Fig. 1. Evolution of the shock structure in the shock frame as a 25 km s" 1 C-type shock interacts with a inhomogeneity of increasing 
density (from ft H = 10 4 cm~ 3 to = 10 5 cm~ 3 ). The left panel shows the neutral density and the right panel the grain-neutral relative 
speed. The profile for the initial C-type shock is given by the thin solid line, while the dashed line shows the shock structure 9.75 
xlO 3 yr after the shock first interacts with the perturbation, the dotted line after 2.56 x 10 4 yr, the dash-dotted line after 6.37 x 10 4 
yr, and the thick solid line the final steady C-type shock structure (reached after 7.92 x 10 4 yr). 



The density and velocity plots in Fig. Q] seem to suggest that 
the neutral flow in sufficiently upstream parts of the precursor 
is quasi-steady. The precursor shows some small oscillations in 
both width and gradient, but these are associated with numerical 
uncertainties. The quasi- steady evolution does not start imme- 
diately after the shock interacts with the density perturbation, 
but only commences after 7 x 10 3 yr. We note that this is roughly 
(L s +L p )/V s with L s the initial shock thickness and L p the pertur- 
bation width. It is the time needed for the shock and perturbation 
to cross each other. 

Shocks with different propagation speeds or moving into re- 
gions with different density contrasts have behaviours similar to 
that of the model shock described above. Thus, such simulations 
give us some idea about the relative strengths of the reflected- 
transmitted shock pairs and the timescales of the J to C-type 
transitions for cores in which the shock moves into a region of 
enhanced density. 

For instance, for a lower density contrast across the pertur- 
bation, the shock does not slow down as much as before. For ex- 
ample, the propagation speed of a shock with an initial speed of 
25 km s" 1 is 21.1 km s" 1 in a region with a density contrast of 3. 
The transmitted shock is only marginally weaker than the initial 
C-type shock with peak temperatures of a few times 10 4 K. The 
timescale for the transmitted shock to become a steady C-type 
shock is slightly shorter than before (6.02 xlO 4 yr or about 1.15 
ion shock crossing times). The post- shock density behind the 
transmitted shock is not much higher than the far downstream 
density leading to a weak reflected shock. The peak temperature 
in the reflected shock is only about 300 K, while it is of the order 
of a few times 10 3 K for a density contrast of 10. As this tem- 
perature is below the threshold for electron cooling to be impor- 
tant, the ion and electron temperatures in the reflected shock are 
similar. Furthermore, the grain-neutral relative velocity, which 
roughly scales with the temperature, in the reflected shock is too 
low for grain sputtering to be important, i.e. \v g — v n \ < 2 km s" 1 
(see discussion in Sect.©. 

We also performed a simulation for an initial shock speed 
of 6 km s" 1 and a density contrast of 10. In this simulation 
the shocks were relatively weak. While strong fast-mode shocks 
considerably heat the incoming gas, weak fast-mode shocks do 



not. For a shock velocity of 6 km s" 1 (or Ma ~ 3), the gas tem- 
peratures are below 10 3 K. The peak temperature in the transmit- 
ted and reflected shock are even less and are about 400 and 80 K 
respectively. Note that the ratio of the temperatures is similar to 
the one for the 25 km s" 1 shock. Because of the low tempera- 
tures in the shock, the grains do not get highly charged and their 
Hall parameter is of the order of unity. The grains, therefore, 
move at a speed different from those of the neutrals and of the 
ions and electrons. As the relative grain-neutral velocity is quite 
small («: 1 km s" 1 ), grain sputtering will not occur. Variations 
in the transmitted- shock structure dissappear after 9.20 xlO 4 yr 
(about one ion shock crossing times) when the shock reaches 
steady state. An interesting difference between this model and 
the strong fast-mode shock models is that the transmitted shock 
width is not much smaller than of the width of the initial C-type 
shock. The ratio of the shock widths is 0.62. This is not surpris- 
ing as the fractional ionisation does not change much across the 
density perturbation due to the low temperatures of the shock. 
Another difference is that, for the weak shock case, the trajec- 
tory in the B y - B z phase space has no spiral node. This is be- 
cause, as the grains are not important charge carriers, the Hall 
conductivity does not modify the upstream conditions. 



3.2. Decreasing density 

For this model, a 25 km s" 1 C-type shock propagates out of a 
dense region with nu = 10 5 cm -3 into a more diffuse one with 
uu - 10 4 cm -3 . Such a model is representative of a molecular 
outflow leaving the dense core surrounding a protostellar object. 

Note that the final steady shock for the case to which Fig. [T] 
corresponds, is similar to the initial shock structure here (see 
Figs. |U[5] and 0). Many features of the shock-clump interaction 
are reversed. As the shock moves into the lower density region, 
the shock speeds up and broadens (see Fig. [4]). Furthermore, the 
spiral By - B z feature does not appear. A rarefaction wave prop- 
agates into the far downstream gas. In this case, the rarefaction 
is a slow-mode wave as the magnetic field decreases across the 
structure. As for multifluid shocks, a multifluid rarefaction wave 
also introduces an ion-neutral drift (see e.g. Fig. [5})). However, 
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Fig. 2. x component of the neutrals (solid), ions and electrons (dashed) and the grains (dotted) velocity in the upstream frame at 
different times for a 25 km s" 1 C-type shock propagating in a higher density region (from nu = 10 4 cm~ 3 to nu = 10 5 cm~ 3 ). (a) shows 
the the initial C-type shock, (b) the the shock structure 9.75 xlO 3 yr after the first interaction of the shock with the inhomogeneity, 
(c) after 2.56 x 10 4 yr and (d) the final steady C-type shock (after 7.92 x 10 4 yr). 

Fig. 3. Ion (dashed), electron (dotted) and neutral temperature (solid) and the average grain charge (dash-dotted) in electron charge 
units for the same shock model as in Fig. [2] at the same times. 
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x(1e16cm) x(1e16cm) 

Fig. 4. Similar to Fig.[U but for a shock interacting with a decreasing density perturbation (from nu = 10 5 cm~ 3 to nu = 10 4 cm~ 3 ). 
The thin solid line is the initial C-type shock, the dashed line the shock structure 5 x 10 3 yr after the shock interacts with the density 
perturbation, the dotted line after 2.40 x 10 4 yr, the dash-dotted line after 5.59 x 10 4 yr and the thick solid line the final steady C-type 
shock (after 9.8 x 10 4 yr). 



the relative velocities between the charged particles and neutrals 
is negligible compared to the velocity of the transmitted shock. 

In the transmitted shock a neutral subshock and a precursor 
are initially present. From the precursor, a C-type shock even- 
tually develops. As in the increasing density, the quasi-steady 
evolution starts after an initial adaptation phase of the order 
1.7 xlO 3 yr (= (L s + L p )/V s ). The transmitted shock reaches 
steady state after 9.80 x 10 4 yr, which roughly corresponds to 
1.2 times the ion flow time through the final shock structure. 

As mentioned before, the transmitted shock speeds up to 
« 28.6 km s" 1 or Ma ~ 13. The shock is thus a strong shock 
and the gas temperatures are high within the shock. (I.e. the ions 
have peak temperatures of 10 5 K.) As the high temperatures are 
due to collisions between the different fluids, the ion-neutral and 
grain-neutral relative velocities within the shock structures also 
remain high (see Fig. [4]). Furthermore, the spatial region with 
high relative velocities is larger than for the case for which re- 
sults are shown in Fig.ffl which suggests that grain sputtering is 
going to be more effective than in the initial C-type shock. 

3.3. Sinusoidal density perturbation 

We also study the interaction of a C-type shock with a clump. 
The shape of the clump is given by a cosine wave with a wave- 
length (A) of 1.46 xlO 17 cm, or more specifically 



n H 



COS 



2n(x — xq) 



A 



- 71 



(nontax + n H ,0\ 

A 2 )' 



with xo the start position of the clump. The maximum density in 
the clump is n H ,max = 2 x 10 5 cm -3 and is 20 times higher than 
the background density (az#,o). Th ese values are repre sentative 
for clump parameters measured by lTafalla et al.l (|2004|) . 

Unfortunately, the evolution of the shock structure cannot 
be described by a combination of the previous two models. The 
behaviour of the shock as it moves up the density gradient is 
similar to the early evolution shown in Fig. [T] The dashed line 
in Fig. [7] and in Fig.[8t) and[9]3 show that a reflected shock and 
a transmitted J-type shock form when propagating up the den- 
sity gradient. The transmitted shock slows and narrows. A spiral 
feature develops in the trajectory in B y - B z phase space. 



The motion back to the low density region does not resemble 
the transition of Sect. 13.21 This is because the transmitted shock 
does not become a steady state C-type shock within a few thou- 
sand years. (The time scale for the shock to adapt to the new 
upstream condition is by itself a few thousand years.) Rather, a 
J-type shock moves down the density gradient. As the shock re- 
acts to the changing upstream conditions, it broadens and speeds 
up. The precursor then develops a neutral subshock just as the C- 
type shock did for the case for which Fig. |4] shows results. This 
means that the shock structure simultaneously has two neutral 
subshocks (see Fig. [St). Interestingly, the neutral and charged 
flows do not have the same post-shock velocities behind the sec- 
ondary subshock. However, this is not surprising as the relative 
velocity of the neutral and charged flow in the original precursor 
is non-zero. 

Although these relativ e velocities are well below the sput- 
tering thresholds derived in lCaselli et al.l (1 19971) . other initial pa- 
rameters may yield relative velocities above this threshold which 
extend over a spatial region larger than the initial shock width. 
The relative velocities steadily decrease as both subshocks even- 
tually weaken. The initial one weakens after a few 10 4 yr. 
Simultaneously, a steady C-type shock forms from the precur- 
sor (see Fig. [7]). The timescale for the steady C-type shock to 
arise is about 1.2 xlO 5 yr or 1.2 times the ion flow time through 
the shock. The final C-type shock is, as expected, identical to the 
initial one. 



4. Summary and discussions 

In this paper we present the first time-dependent simulations of 
oblique C-type shocks in dusty plasmas interacting with density 
perturbations. We studied the transient evolution of the shock 
structure for each of several different types of density perturba- 
tions, i.e. increasing, decreasing and clump-like, and examined 
the dependence of this interaction on shock speed and density 
contrast. 

When a steady C-type shock encounters a density perturba- 
tion, the shock reacts to the changing upstream condition by 
breaking up into multiple waves. While a J-type shock moves 
into or out of a dense region, a shock or rarefaction wave, respec- 
tively, is reflected and propagates in the reverse direction. This 
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Fig. 5. Similar to Fig. [2 but for a shock interacting with a decreasing density perturbation (from nn = 10 5 cm~ 3 to nu = 10 4 cm -3 ). 
(a) shows the initial condition, (b) the shock structure 5 x 10 3 yr after the interaction with the density perturbation, (c) after 2.40 x 10 4 
yr and (d) the final steady state shock (after 9.8 x 10 4 yr). 



Fig. 6. Similar to Fig.O but for the shock and times shown in Fig. [5] 
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Fig. 7. Similar to Fig.[T] but for a shock interacting with a clump that has a maximum density of nu = 2 x 10 5 cm~ 3 . The solid line 
shows the initial and final steady shock (after 1.2 xlO 5 yr), the dashed line the shock structure 3.9 xlO 3 yr after the shock interacts 
with the clump, the dotted line after 1.98 xlO 4 yr and the dash-dotted line after 1.46 xlO 5 yr. 



reflected shock or wave is present because the post-shock flow 
behind the transmitted shock is different from the post-shock 
flow of the initial shock. The relative strength of the transmit- 
ted and reflected components depends on the initial shock speed 
and the density contrast of the perturbation, as these control 
the post-shock properties of the transmitted shock, e.g. a lower 
density contrasts leads to a smaller difference between the far 
downstream flow and the transmitted post-shock flow. Hence, a 
weaker rarefaction wave or shock is needed to adjust this change. 
For the interaction with a clump, where the density decreases 
again, the situation gets more complicated as the transmitted J- 
type shock forms a second sub shock in the neutral flow of the 
precursor. It is thus important to realise that we expect a range of 
different J- and C-type shocks and multifluid rarefaction waves 
in inhomogeneous clouds. 

In our models, we have focused on the evolution of the trans- 
mitted shocks. Such a shock is initially J-type and contains a sub- 
shock in the neutral flow. It again becomes a steady C-type shock 
but with a different shock width and velocity. The timescale for 
it to become steady is of the same order as the ion-flow time 
through the final shock structure, i.e. about 0.5-1.5 xlO 5 yr. An 
important consequence of these timescales is that, if a shock in- 
teracts with density inhomogeneities on timescales shorter than 
this, no steady C-type shock exists. All shocks are thus J-type, 
although some shocks will have weak subshocks and are then 
only marginally distinguishable from C-type shocks. 

Our simulations show that the steady dusty C-type shock 
develops from a shock in which the sufficiently far up- 
stream parts of the precursor are steady. Such an evolution 
was also obser ved for multifluid shocks in dustl ess, weakly - 
ioni sed gases ([Flower & Pineau des Fqretsl 1 19991). Following 
the iFlower & Pineau des Foretsl (1 19991) result. ILesaffre et alJ 
(120041) " developed a quasi- steady method to follow the temporal 
evolution of multifluid shocks. This method involves the treat- 
ment of an evolving non- stationary J-type shock as a sequence 
of truncated C-type shocks each of which has a neutral subshock 
at the point where the ion-flow time corresponds to the age of 
the shock. While this approach was only validated fo r shocks in 
dustless, weakly ionised gases (ILesaffre et all 120041) . our result 
confirms quasi- steady models can also be used for dusty shocks 
and justifies the use o f such method s for t he temporal evolu- 
tion of dusty shocks by lGusdorf et al] (120081) . However, we also 



find some limitations of the quasi- steady approach when applied 
to the interaction of shocks with density perturbations. Firstly, 
this method can only be used after an early adaptation phase in 
which the shock adjusts to the changing upstream conditions. 
Secondly, this approach is only valid for a C-type shock inter- 
acting with a density perturbation. For a J-type shock, a second 
neutral subshock forms in the precursor. Behind this secondary 
shock, the velocities of the neutral and charged fluids are not 
the same. Such a situation, which is likely to occur in a clumpy 
medium, cannot be modelled with a quasi- steady approach. 

In quiescent cold clouds, silicon is nearly depleted from the 
gas phase and stored in the core of dust grains. However, ob- 
servations of SiO emission near protostellar objects show that 
some silicon is returned to the gas-phase. It is thought that this is 
due to s puttering of the dust grains a nd grain-grain collisions in 
shocks dMartin-Pintado et all Il992l) . Although these processes 
are not included in the current code, the change of the grain- 
neutral relative velocity serves as an indication. Grain mantle 
sputtering requires a grain-neutral relative velocity of « 6 km s" 1 
whereas core sp uttering requires a m uch higher relative velocity 
of « 19 km s" 1 dCaselli et all 1 19971) . For our model parameters, 
with a shock speed of 25 km s" 1 , the threshold for core sput- 
tering is only just attainable. Therefore, we do not expect much 
SiO emission from these models. Howeve r, grain mantle sput- 
tering occurs in all the initial C-type shocks. DTmenez-Serra et all 
( 2008) have suggested that if the grain mantle consists of a small 
fraction of silicon, sputtering by heavy molecules such as CO 
could erode the mantles even in low velocity shocks and account 
for the narrow SiO lin e emission observed in som e young out- 
flows (e.g. L1448-mm Jimenez- Serr a et all 120041) . Silicon ero- 
sion from the mantle would be saturated for shock velocities be- 
tween 10 and 20 km s" 1 . This implies that the SiO enhancement 
from mantle erosion is saturated in all our models. Variations of 
the SiO emission then arise as the spatial area over which sput- 
tering occurs changes (i.e. the shock width increases/decreases). 
Furthermore, we expect a small contribution to mantle sputter- 
ing from reflected waves. In subsequent papers we will include 
grain mantle and core sputtering as well as grain-grain collision 
terms to describe the release of SiO from the gra ins to the gas 
phase and use this to calculate t he SiO emission dCaselli et all 
119971: iJimenez-Serra et al.l.l2008l) . 
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Fig. 8. Similar to Fig. [2 but for a shock interacting with a clump that has a maximum density of nu = 2 x 10 5 cm 3 . (a) shows the 
initial and final shock structure (after 1.2 xlO 5 yr), (b) the shock structure 3.9 xlO 3 yr after the shock first interacted with the clump, 
(c) after 1.98 xlO 4 yr and (d) after 1.46 xlO 5 yr. 

Fig. 9. Similar to Fig. [2 but for the shock and times shown in Fig. [SI 
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